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Abstract 

A Computational Fluid Dynamics method is introduced to study the nucleation, coagulation, evaporation and condensation of poly¬ 
disperse particles in multi-dimensional multi-species laminar and turbulent flows with strong mass and energy coupling between the 
phases. The model is based on the Reduced Navier-Stokes (RNS) methodology and incorporates a lognormal aerosol moment method 
to describe the evolution of suspended particulates. The model is validated against available analytic and numerical techniques for one 
and two-dimensional geometries. The sensitivity of both the bulk flow and aerosol particle properties to boundary layer effects, even in 
high Reynolds number flows, is demonstrated for transonic two-phase flow in wet steam de Laval nozzles. It is shown that the developed 
model and method are useful for the prediction of particle properties such as mass and number concentrations, geometric mean particle 
size and standard deviation of the particle size distribution in multi-dimensional turbulent compressible aerosol flows with strong heat 
and mass transfer. 

© 2006 Elsevier Ltd. All rights reserved. 


1. Introduction 

Understanding and predicting polydisperse aerosol 
behavior is critical to the development of many advanced 
technologies related to combustion, pollution prediction 
and control, fouling, powder and materials manufacture 
and even to the administration of drugs and purification 
of silicon chips. In many of these systems, the bulk fluid 
behavior and that of the suspended particulates can 
become strongly coupled due to mass, momentum and 
energy transfer between the aerosol vapor and condensed 
phases. Understanding such multi-phase processes requires 
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detailed knowledge of complex fluid behavior coupled with 
transport and dynamics of gas species and liquid and solid 
particulates. Aerosol parameters such as size, polydisper- 
sity, concentration, composition and morphology are a 
function of residence time in critical regions of the flow, 
local thermodynamic conditions such as temperature and 
pressure (and hence saturation conditions), as well as mix¬ 
ing efficiencies and chemical reactions affecting heat, mass 
and momentum transfer. These are in turn functions of fac¬ 
tors such as entrainment of surrounding fluid due to viscos¬ 
ity and large and small scale structures in the flow such as 
vorticity, recirculation and turbulence. In many of these 
circumstances, the coupling between the gas phase and sus¬ 
pended aerosol particles can become significant and 
neglecting such effects can lead to erroneous predictions 
of overall system behavior. 



D. P. Brown et al. / Computers & Fluids 35 (2006) 762-780 


763 
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kih. volume moment of the distribution (m 3k /m 3 ) 
mixture molecular weight (Ref. [11]) 

(Ef=i PsCs + P a Ca)/Ef=i (kg/mol) 
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y mixture ratio of specific heats (1 — R/MWC p ) 

(ND) 

yk kih coagulation coefficient (ND) 

3 boundary layer thickness (m) 

A change due to aerosol formation and growth 
(ND) 

s k Lennard-Jones Temperature (K) 

k thermal conductivity (kg/s 3 K) 

rj collection Efficiency (ND) 

X mean free path of gas mixture! X* = 

(m) f B 
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<P energy dissipation function (kg/m s) 

p viscosity (kg/m s) 
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space (ND) 
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Subscripts 

a aerosol phase 

j €,ri,t 

Con due to condensation 

C convective 

D diffusional 

f bulk fluid 

g geometric mean 

I inertial 

k Mi aerosol moment 

1 laminar 

L local 

Nuc due to nucleation 

p particle 

Rxn due to chemical reaction 
s gas phase species s 

t turbulent, total 

T thermophoretic 

v viscous 

w wall value 

Z surface tension 

rj, £ differentiated with respect to £, rj , £ 

1 monomer 

32 Sauder mean 

oo initial, free stream or characteristic value 

Superscripts 

* dimensional quantity 

0 critical value 

i x,y,z 

c continuum 

fm free molecular 

k to the power k 


The computation of such complex processes can be car¬ 
ried out using a variety of numerical techniques but has 
typically been approached by solving the flow in an Eule- 
rian reference frame and tracking particles through the 
flow field in a Lagrangian reference frame. In such 
schemes, aggregate aerosol properties are computed by 
summing contributions of discrete particle sizes and com¬ 
position from the complete size distribution within a given 
flow or by integrating particle properties along streamlines. 
However, in order to accurately account for the particle 
behavior (e.g. agglomeration, nucleation, and growth) 
and thus its influence on the bulk fluid, the population bal¬ 
ance equation (PBE) [1-3] should be solved together with 
the continuity equation, momentum balance equations 
and the energy balance equation. To date, attempts to 
incorporate the PBE into computational fluid dynamics 
(CFD) codes have been limited. A logical extension of 
the Eulerian-Lagrangian approach is the Monte-Carlo 
simulation technique in which the solution of the PBE is 
formulated in terms of its stochastic equivalent [4,5]. Here, 


a representative population of particles evolves according 
to an appropriate balance of probabilities. However, due 
to the large number of particles that need to be tracked 
in order to provide a statistically reasonable representation 
of the physics, the computational expense is still too large 
for the formulation to be implemented in practical 
applications. 

Significant reductions in computational effort can be 
achieved using an Eulerian formulation for multi-phase 
processes [6] and, in recent years, several researchers have 
tried to couple the PBE via Euler representations of the dis¬ 
persed phase. In general, these can be classified as mono- 
disperse methods, discrete population balance (DPB) or 
sectional methods and moment methods. In monodisperse 
methods, details of the PBE are neglected and only the 
overall volume fraction of the dispersed phase is accounted 
for [7-9]. In the DPB technique [10], the continuous parti¬ 
cle size distribution is represented as a collection of discrete 
particle size bins each representing a range of particle prop¬ 
erties and with appropriate transport and dynamics terms 
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applied for each size class. There have been numerous 
models developed for one-dimensional flows using such 
techniques e.g. [11-18]. In order to achieve reasonably 
accurate resolution and to overcome numerical diffusion, 
however, the particle size distribution must be divided into 
a large number of bins. Consequently, implementation in 
CFD codes becomes prohibitive for all but the most restric¬ 
tive cases [19-26]. 

An alternative approach is the methods of moments 
(MOM) in which the moments of the particle size distribu¬ 
tion are solved by integrating out the internal coordinate 
(e.g. diameter or volume). This allows the problem to be 
reconstructed in a reduced variable space and in which 
each moment is a transportable scalar, however, the prob¬ 
lem of mathematical closure arises [27,28]. Three methods 
are typically employed to resolve the closure issue: 
through the assumption of a closed form of the distribu¬ 
tion function in which the properties of the distribution 
are a function of the solved moments [29-31], though 
complex closure rules based on interpolation between 
moments [28] or through quadrature methods [32]. 
Depending on the method of closure, the number of 
unknowns can be reduced to as few as three, which makes 
implementation into CFD codes tractable [33-39]. How¬ 
ever, to date, many issues related to interphase coupling 
have been ignored. In particular, the combined effects of 
mass and heat transfer with the evolution of the polydis- 
perse aerosol size distribution has received no attention 
in the literature. 

Our previous work has modeled polydisperse aerosol 
transport in two and three-dimensions as well as the effects 
of chemical reaction and gas phase species mixing on the 
evolution of the aerosol size distribution [40-42]. The pres¬ 
ent work describes a method of computing aerosol systems 
with strongly coupled heat and mass transfer and extends 
many aspects of previous research into a comprehensive 
model of coupled fluid/species/aerosol behavior which 
can be applied to a wide range of problems. This method 
uses an Eulerian moment form of the general dynamic 
equation for aerosol behavior developed using a lognormal 
size distribution in conjunction with a Reduced Navier- 
Stokes (RNS) formulation [43] for bulk fluid modeling 
which allows large time steps and rapid convergence. It is 
demonstrated that this RNS/moment model is effective in 
solving for behavior of polydisperse aerosols in general 
laminar/turbulent, incompressible/transonic/supersonic, 
flows with strong coupling between the bulk flow, individ¬ 
ual gas phase species and suspended particles. Many 
important parameters, such as spatial/temporal distribu¬ 
tions of aerosol properties and depositional efficiencies, 
can be calculated directly without having to sum up results 
for monodisperse aerosols over the entire particle spectrum 
or without having to solve for a large number of 
unknowns. This formulation can significantly reduce com¬ 
putational effort as compared with schemes based on 
Lagrangian or Eulerian sectional representations of aerosol 
behavior. 


2. Governing equations and boundary conditions 

Four sets of governing equations are considered here: 
those related to the transport of mass, momentum and 
energy in bulk fluid flow, those related to the conservation 
of individual species in the fluid, those related to the pro¬ 
duction, dissipation and transport of turbulence, and 
those related to the transport and dynamics of suspended 
particles within the fluid flow. These sets of equations can 
be strongly coupled to each other and must be solved 
simultaneously in many instances, though, in the low par¬ 
ticle concentration limit, the coupling between the fluid 
dynamics and particle dynamics is one-way. Details of 
the derivations of the governing equations can be found 
in [40-42]. 

2.7. Bulk fluid equations 


The bulk fluid flow is described by the compressible 
Reduced Navier-Stokes (RNS) [43] equations. They are 
derived using a perturbation expansion for large Reynolds 
number: omitting higher order stream wise diffusion terms 
in the full Navier-Stokes equations. The formulation 
reduces to the Euler equations when all viscous diffusion 
and heat transfer terms are neglected. This results in a 
composite system of the full Euler and higher order 
boundary layer equations. The RNS equations can thus 
be used to solve both the inviscid and viscous regions 
simultaneously, providing an efficient way of solving both 
internal and external flow fields at a lower computational 
cost than many full Navier-Stokes (FNS) solvers. Other 
advantages of the RNS method include simplified 
consistent boundary conditions, sharp shock capturing 
without the use of explicit artificial viscosity, and faster 
convergence than time-dependent Navier-Stokes solvers. 
The Full Navier-Stokes solution can be recovered by 
replacing the remaining diffusion terms in a deferred- 
corrector (DC) [44]. 

The RNS method has been verified extensively over a 
wide range of conditions including 1-D/2-D/3-D, incom¬ 
pressible/hypersonic, viscous/in viscid, laminar/turbulent, 
steady/unsteady and internal/external flows. The present 
code is written in conserved variable form and incorporates 
velocity flux splitting for continuity and convective terms. 
Using the following non-dimensionalization, 
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(1) 


the RNS equations with source terms due to fluid/aerosol 
particle coupling and with multiple species in the gas phase 
can be written in vector form as 
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where p is the mixture density where (p = p f + p' a ) and 
(p f = Y^LiPs)') diag * s a diagonal matrix; t/ f , F f and Up 
are contravariant velocities. w f , r f and w f are Cartesian 
velocities where 


Wf = XgUf +Xr,Vf (3) 

Vf =y^Uf +y n v f + 

Wf = z^C/f + z^Cf + z^Wf 


H t is the total enthalpy. A H t and A m are the heat and mass 
addition due to aerosol formation and growth, qj the heat 
flux. Sources due to momentum transfer and buoyancy are 
not considered here, x, y and z are Cartesian coordinates. 

r\ and C are transformed coordinates, g is the inverse 
Jacobian of the coordinate transformation, t-, <Pj represent 
the RNS approximated viscous stress tensor and viscous 
energy dissipation respectively. Subscripts connote the 
derivative of the subscripted quantity with respect to 
the subscripted variable, y is the ratio of specific heats of 
the aerosol and M is the Mach number for the mixture de¬ 
fined from the relations [11] 


C p — C s C Ps + c a c, a 


MW= + 
7 = (1 - tf/MWC,,) 


Ns n 

P s Cs 
^ MW, 

s= 1 A 


(4) 

(5) 

( 6 ) 


M = 



+ V * 2 + w* 2 
yRT* 


(7) 


where N s is the number of gas phase species, C p is the spe¬ 
cific heat at constant pressure and MW is the molecular 
weight for the mixture, including the condensed phase. p a 
and C Pa are calculated based on mass concentration of aer¬ 
osol per unit gas volume. Note that equivalent molecular 
weight for the aerosol, MW a , is neglected in the MW calcu¬ 
lation since it appears in the denominator of a summation 
and is much larger than the gas phase molecular weights. 


Closure of the equations is achieved with the equation of 
state 

7co Ml2 a P = p (h x - 1 -( 7oo - 1 )M 2 Ju 2 + v 2 + w 2 )^ (8) 

2.2. Aerosol equations 

The present method uses an Eulerian moment form of 
the general dynamic equation (GDE) developed assuming 
a lognormal size distribution [29,45] which has been found 
to exist in a wide range of aerosol systems [46]. The model 
incorporates contributions for particle transport (due to 
convection, diffusion, thermophoresis and inertia) and 
dynamics (due to coagulation, nucleation, evaporation 
and condensation). 


2.2.7. Aerosol transport 

The aerosol conservation equation can be written as 

+ v . (« p (i; p )Up) = 0 ( 9 ) 

where n p (v p ) is the dimensionless particle size distribution 
function, U p is the dimensionless particle velocity vector. 
Up can be defined as U p = U c + U D + Uj + U T where 
U c = convective velocity, U D = diffusion velocity, 
Ui = inertial velocity, and U T = thermophoretic velocity. 
An overview is given here. Convective velocity is given by 
the flow field. Diffusional, inertial and thermophoretic 
velocities are strong functions of the flowfield characteris¬ 
tics such as velocity gradients, concentration gradients 
and temperature gradients as will as local thermodynamic 
conditions. 

The diffusional velocity is given as U D = —Z>*(Vln(w*)) 
where D* (the particle diffusion coefficient) is derived from 
kinetic theory as D* = k B T^T*/m p [47]. t* is the particle 
relaxation time for particle motion due to the effects of 
the viscous drag force. For sub-critical Stokes number 
(Re p < 1), t* is given by 



where C c is the form of the Cunningham correction factor 
given in [45]. Inertial velocity is a modification of the form 
given by [49], Ui = — r*(UJ • V)U£. The steady state conser¬ 
vation equations can be written as 


(u;-v)u; + V;-u*) + 2v- J p; = o 

T P p 

v • (p;u;) = o 
p> rj wn 

p m* 
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( 10 ) 


where t is the relaxation time appearing the definition of 
the particle drag 
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The current scheme assumes a relation first order in t, thus 
U p =U f + U D +U I + U T (12) 

where 

Ui«(U f -V)Uf (13) 

Several forms of the thermophoretic velocity have been 
given in the literature [47,48,50,45]. In general they can be 
written in the form U j =—K T m v \7T* { / where K v = 
K T (/P n Kn p ) and Kn p is the Knudsen number of the particle. 
For the present cases, K T is assumed to be a linear function 
of since it is generally a weak function of d*, and 2*. 
Using the following non-dimensionalization in conjunction 
with the non-dimensionalization used for fluid variables 


Ml 
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V P — 4 3 

- nr- 5 
3 IU poo 


and 
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the particle velocity in non-dimensional form is 
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where Re is the fluid Reynolds number 
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(14) 


(15) 


The Peclet number is defined as 


Pe = ScRe L = 


D* 


Pf Pp UqqLqO 

m Pf 


2.2.2. Aerosol dynamics 

The general dynamic equation (GDE) including Brown¬ 
ian coagulation, condensation, evaporation, and nucle- 
ation is given by [50] 


p»; 8 (gq 

D t 0r* 




X [ P f 2 (P,vl,%,n*,r)dv*dv* 

Jo 


(16) 


where the first term on the left is the total rate of change 
(substantive derivative) in number concentration including 
the non-convective transport terms, the second is the rate 
of change in number concentration due to condensation 


at rate G, the third is the rate of formation of new particles 
of critical volume v® at rate / Nuc . The RHS represents the 
coagulation rate, where /? is the collision frequency func¬ 
tion. From classical nucleation theory [50], the nucleation 
rate for a given gas phase species (s) is given by 


r 2 (h n 

7 Nuc s — F^ ucs n s S\ s I 

V 




Z; exp 




where 


pO _ n ( 4 K Z Z S \ 3 
K s ~ 6 \\n(S s )) ' 


(17) 


Since there is significant uncertainty in the prediction of 
nucleation rates, it is common to tune the classical rate to 
match data. Two methods are typically employed: a pre¬ 
multiplication factor for the entire nucleation term or an 
adjustment to the surface tension [51]. In the present work, 
the nucleation rate has been scaled by ^nuc an d surface 
tension is calculated based on a fourth order curve fit to 
experimental data. is a parameter that adjusts for the 
ratio between the surface tension on a flat surface and that 
of a nucleating droplet. 


2.2.3. Polydispersity modeling 

There exist a number of methods to describe the behav¬ 
ior of polydisperse aerosols. The most general method 
assumes nothing about the aerosol distribution and models 
the general dynamic equation (GDE) using direct numeri¬ 
cal simulation of individual particles. The computational 
expense of such simulations, however, means they are gen¬ 
erally beyond the capabilities of current supercomputers, 
even for rudimentary problems and alternate methods 
must be used for practical applications. A typical approxi¬ 
mation employed is the discrete sectional method [10] 
which sacrifices some accuracy in favor of computational 
efficiency. In this method, the size distribution is divided 
into a finite number of bins in which properties are aver¬ 
aged. The accuracy of the method is a function of the res¬ 
olution of the particle size bins. Polydisperse aerosol 
behavior can be modeled more efficiently by making an 
assumption about the shape size distribution and solving 
the GDE in moment form for the properties of that distri¬ 
bution. This is the approach taken here. 

It is assumed that the particle size distribution is log¬ 
normal, i.e. 

«p( u p) = 1/(3\/27u>* In< t) exp{-ln 2 (D;/Dp/181n 2 (r} (18) 

where r* and r* are particle volume and geometric mean 
particle volumes, and is the geometric standard deviation 
of the size distribution. The lognormal distribution has 
the following properties: 


v*7 = vf exp (^/v 2 ln 2 <r) 
M*=M>fex pjJ^lnV) 
where M* k as J 0 °° v* k n* p (v* p , f) cb*. 


( 19 ) 
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In general, three equations are required to define the 
lognormal distribution. These can be obtained by multiply¬ 
ing the conservation equation by v^ k (where k is the 
moment of the size distribution) and integrating over all 
aerosol volumes to give the moment form of the particle 
conservation equation. We have chosen to solve for the zer¬ 
oth, first and second moment conservation equations by 
setting k — 0,1,2. Other choices for the solved moments 
are possible. The geometric mean particle volume (r g ) 
and the standard deviation of the distribution (cr) can be 
defined from the zeroth, first and second moments as 


Vo = 


M\ 


M l 0 5 M° 2 5 ' 


Inn = - In 


MqM' 

m\ 


ln 2 c 


( 20 ) 


where is the reference geometric mean standard 
deviation of the distribution [45]. 

In general, G and /? take different forms for continuum, 
slip, and free molecular regimes [55,56]. One set of equations 
can be used for all regimes by using a harmonic average of 
the two equations for continuum and free molecular systems 
[29]. Substitution of U p into the GDE yields a moment form 
representing three simultaneous partial differential equa¬ 
tions. The properties of the distribution can be used to 
arrange the equations such that only the ki\\ moment 
appears explicitly only in the kih moment equation [42], thus 
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( 21 ) 


where k = 0,1,2. This resulting system of equations is 
strongly coupled both to each other and to properties of 
the flow (through Sc k , y k , oc k , I s , 7>, U f , C s and S s ). 

The present model incorporates the factor ^cond to 
adjust for the rate at which enthalpy of phase change can 
be conducted to the surrounding fluid. This is analogous 
to the condensation coefficient ot c used by [8] and the Knud- 
sen correction </> used by [34] for the continuum regime. 
The polydisperse Stokes and Schmidt numbers (St k and 
Sc k ) can be written as 

Sc 0 = \/{A 3 v~ x ^ 3 exp(0.51n 2 cr) -\- A 4 v~ 2 ^ 3 exp(2ln 2 cr)) 

Sci = 1/(A 3 v ~ 1/3 exp(-2.51n 2 cr) + A 4 v~ 2 / 3 exp(-41n 2 a)) 

Sc 2 = l/(A 3 v ~ 1/3 exp(-5.51n 2 <x) +A 4 v~ 2 / 3 exp(-101n 2 <x)) 

( 22 ) 

St 0 = A\v 2 ^ 3 exp(21n 2 a) + A 2 v l J 3 exp(0.51n 2 a) 

Sti = A\v 2 J 3 exp(81n 2 cr) + A 2 v l J 3 exp(3.51n 2 a) (23) 

St 2 =Aiv 2 J 3 exp(141n 2 o) + A 2 v l J 3 exp(6.51n 2 cr) 
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Coagulation and condensation coefficients y K and a K are 
calculated from the harmonic average of the limiting cases 
for small (free molecular) and large (continuum) aerosols. 
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(26) 


The constants B u B 2 , B 3 , B 4 and B 5 are related to thermo¬ 
dynamic properties of the system [29]. b Q and b 2 are given in 
[55]. 


2.3. Species conservation equations 

For multi-species applications, the transport and con¬ 
sumption of individual species must be considered. In gen¬ 
eral, a separate continuity equation is required for each gas 
phase species. Thus for all species s, the saturation equa¬ 
tion takes the form 


6pC, 

dt 


+ V- 
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where C s is the gas phase mass fraction for species s. The 
RHS represents the sources in mass due to nucleation, 
evaporation, condensation and chemical reaction. 
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2.4. Turbulence model equations 


For turbulent calculations, a Reynolds averaged form of 
the flow, species and aerosol equations is used where the 
turbulent viscosity is calculated from a turbulence model 
from the Baldwin-Barth turbulence model. The resulting 
turbulent viscosity is then used to calculate the turbulent 
mass diffusivity in the aerosol moment and gas phase spe¬ 
cies equations as well as the total viscosity in the fluid equa¬ 
tions. The Baldwin-Barth [52] one equation turbulence 
model is suitable for both internal and external flows. 


9(vi Ret) 
3 1 


■ U • V(vRe t ) 
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2.4.1. Turbulent contribution to aerosol moment equations 

The state of the art for turbulent modeling of particulates 
is weak, though several attempts at studying the effects have 
appeared in the literature. In general, it has been observed 
that particles in the micron and submicron range experience 
turbulent fluctuation of the same order of magnitude as the 
bulk fluid [53]. Larger particles do not experience a signifi¬ 
cant amount of turbulent fluctuation. 

Reynolds’ analogy between heat and mass transfer 
assumes D t /oc t = 1 [54]. A similar analogy can be written 
for viscous and mass diffusivity, i.e.: 0 <v t /D t < 1 which 
can be written as a function of the particle Stokes number 
where here the turbulent flow eddy diffusivity is calculated 
using the Baldwin-Barth Turbulence model. The effects of 
the turbulent diffusivity are then added to the diffusion 
term via a correction factor to the viscosity. 

Turbulent Correction Factor = 1 + — (29) 

P\ W 

where ^ = f(St k ). When Stokes number is large, j[St k ) —► 
0. When Stokes number is small, f[St k ) —> 1. In this work, 
particles are in the submicron range and it is safe to assume 


a small Stokes number as they are transported effectively 
with the turbulent eddies. 


2.4.2. Coupling source terms 

Gas and aerosol phases are indirectly coupled through 
thermodynamic parameters such as pressure and tempera¬ 
ture, as well as directly influenced through source terms 
due to mass and energy transfer between the phases. Cou¬ 
pling source terms for mass and energy transfer between 
the bulk flow, individual gas phase species and aerosol 
moments are determined from the relevant system para¬ 
meters via 
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(30) 


where Am and AH t appear in the bulk flow continuity and 
energy equations respectively and A m s appears in each 
species equation. The first term on the RHS represents 
the mass exchange due to nucleation and the second term 
represents the change due to condensation/evaporation. 


2.4.3. Boundary conditions 

The boundary conditions for the bulk flow are: no-slip 
at walls (pressure calculated) and P = constant at free sur¬ 
faces (normal velocity calculated); H t specified at the inflow 
and boundaries; inflow density and velocity are specified; 
P = constant at the inflow for supersonic inflows. Pressure 
alone is specified at the outflow for subsonic flows and d 2 P/ 
6£ 2 = 0 for supersonic outflows. For subsonic inflow, a 
pressure and density are calculated based the characteris¬ 
tics of the Navier-Stokes equations via 


p=pjr~^ l) 

p = Pt/f 


(31) 


where 

C + = 1 + 2/(y — 1)^00 
C~ = Uf + 2/ (7 — 1 )a 

f + - = 1 +1(7-1) 

Boundary conditions for the aerosols are; M 0 , and 
M 2 = constant at the inflow, dM k /drj = 0, dM k /d{, = 0 (slip 
or zero flux) or M k = constant (non-slip) at the walls, and 
zero flux at free surfaces. The value of the constant can be 
used to reflect the nature of the boundary. M k = 0 reflects 
total capture of all particles that reach the wall. M k 
between zero and unity reflects partial capture of particles. 
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Boundary conditions for the gas phase species are: 
C s = constant at the inflow and dCjdrj = 0, dCjd£ = 0 
(zero flux) at the walls and free surfaces. Boundary condi¬ 
tions for the turbulent viscosity are: vRt = constant at the 
inflow and free boundaries and vRt = 0 at the walls. 

In all equations, for reversed flow at the downstream 
boundary, the flare approximation (in which streamwise 
convective terms are neglected) is used to eliminate down¬ 
stream velocity, moment, species and turbulence boundary 
condition requirements. For 1-D and 2-D cases, appropri¬ 
ate cross-flow derivatives are dropped. 

3. Discretization and solution procedure 

Flow, species, turbulence and aerosol behavior are com¬ 
puted in separate modules using a finite difference discreti¬ 
zation scheme and simultaneously marched in space based 
on the principles of the RNS methodology [43]. Details of 
the discretization of the governing equations and the local 
and global solution procedure are described below. 


The equations are centered in the /-plane and staggered 
in the cross-plane (j,k- plane). Figs. 1 and 2 show the cross 
plane equation locations and grouping for internal and 
external flows. Note that the selection of variables associ¬ 
ated with each equation is different than used in the tradi¬ 
tional RNS methodology. In general, on solid walls, the 
normal momentum equation is used to calculate pressure 
and the normal velocity is prescribed. For free surfaces, 
the normal velocity is calculated and the pressure is 
prescribed. 

For two-dimensional cases, the stencils reduce to three 
planes in the ^-direction (Fig. 3). Boundary conditions on 
the k— 1 and k = 3 faces are slip conditions (0/0£ = 0). 
The 3 xy max cross plane is then inverted as in the full three 
dimensional problem. For axisymmetric cases, the three 
planes are placed at an angle with a common vertex at 
y'=l and a small radius is added at the centerline to avoid 
a numerical singularity when r = 0. 

3.2. Global and local solution procedure 


3.1. Discretization 

For all sets of equations, parabolic convective terms are 
upwind differenced with respect to the dominant flow direc¬ 
tion (£), and central differenced in the cross flow directions 
(f/,Q. Elliptic diffusion, heat flux and dissipation terms are 
central differenced in the cross flow directions and neglected 
in the streamwise direction. Additional terms that appear in 
the aerosol moment equations, i.e. thermophoretic and 
inertial terms are also treated according to their character¬ 
istics. Elliptic thermophoretic terms are neglected in the 
axial direction and central differenced in the cross plane. 
Parabolic inertial terms are upwind differenced in the axial 
direction and central differenced in the cross plane. 

Since it is desired to determine the flow field in arbitrary 
geometries, and since it is necessary to resolve the solution 
in zones of high gradients, the finite difference representa¬ 
tions of the governing equations are written in generalized 
three-dimensional non-orthogonal curvilinear coordinates 
and are solved in delta form to reduce the number of 
matrix inversions. In three dimensions the equations are 
discretized according to Table 1. 


In the RNS methodology, axial diffusion terms are 
neglected, which allows an efficient spatial marching proce¬ 
dure in which only a single axial plane is treated implicitly. 
Coupling effects are calculated during the local solution 
procedure. Globally, a multi-sweep, space marching (global 
pressure relaxation) procedure is used in which the solution 
is marched in the approximate flow direction (£ or /-direc¬ 
tion). The full set of equations is converged at each £ station 
before moving on to the next. Each /-plane is solved implic¬ 
itly, / + 1 and / - 1 planes are treated explicitly and globally 
relaxed. Pressure is flux split in the streamwise direction 
according to the streamwise component of Mach number 
based on the parameter co = min(M^, 1) where M includes 
the particle contribution. 
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Table 1 

Equation differencing scheme 


Equation 

Equation location 

Differencing 



/-direction 

y-direction 

/^-direction 

Continuity 

1 

<N 

1 

Upwind 

Trapezoidal 

Trapezoidal 

“<T Momentum 

(Uj,k) 

Upwind 

Central 

Central 

“rj” Momentum 

{i,j+ 1/2, fc) 

Upwind 

Trapezoidal 

Central 

Momentum 

(i,j,k+ 1/2) 

Upwind 

Central 

Trapezoidal 

Energy 

O'J.fc) 

Upwind 

Central 

Central 

State 


NA 

NA 

NA 

Aerosol moments 

( i,j,k ) 

Upwind 

Central 

Central 

Gas phase species 

(i,j,k) 

Upwind 

Central 

Central 

Turbulence model 

( Uj,k ) 

Upwind 

Central 

Central 
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Fig. 1. Cross-plane (jk) equation locations and groupings for 3-D internal 
flows. C = continuity, £ = ^-momentum, r\ = ^-momentum, £ = (-momen¬ 
tum, be = boundary condition. Subscripts denote the unknown associated 
with each equation in coefficient matrix. Fluid state and energy equations 
as well as all moment and species equations are located at the grid points. 



Fig. 2. Cross-plane (jk ) equation locations and groupings for 3-D internal 
flows. C = continuity, £ = ((-momentum, rj = ^-momentum, £ = (-momen¬ 
tum, be = boundary condition. Subscripts denote the unknown associated 
with each equation in coefficient matrix. Fluid state and energy equations 
as well as all moment and species equations are located at the grid point. 


Streamwise convective terms in all the equations as well 
as inertial terms in the moment equations are upwind dif¬ 
ferenced based on the sign of t/ f . In unseparated flow, 
upstream elliptic influences appear only through the pres¬ 
sure in the RNS equations (and, consequently, through 
the saturation ratio in the species and aerosol equations). 
Thus multiple sweeps are performed to carry the down¬ 
stream pressure information forward. In reversed flow 
regions, all upwinded variables are relaxed in a similar 
manner since velocity at i + 1 must then be used in the 
/-plane calculations. For fully supersonic flows without 


recirculation, the governing equations become parabolized 
and only a single global sweep is required to achieve 
convergence. 

Coupling between the modules occurs during the local 
solution procedure. At every /-station, the sets of governing 
equations are solved sequentially. Each of these equations 
is nonlinear and thus must be iterated to convergence using 
Newton’s method. Additionally, since flow, species, turbu¬ 
lence and aerosol equations (including local particle veloc¬ 
ities) are also coupled to each other in a non-linear fashion, 
a second level of iteration is required. 

In the RNS equations and species equations, the com¬ 
plete system is inverted at once. The three moment equa¬ 
tions are solved in a semi-coupled fashion, in which the 
three equations are considered simultaneously with the 
coefficients on the coupling terms lagged one iteration. It 
has been found that this increases the stability over solving 
the equation simultaneously since the source terms in the 
moment equations due to nucleation, evaporation, conden¬ 
sation and coagulation can vary over many orders of mag¬ 
nitude. Increased stability is achieved at the expense of a 
slightly reduced rate of convergence. 

3.3. Matrix inversion routines 

The discretized fluid equations results in a block penta- 
diagonal coefficient matrix which is inverted using a sparse 
matrix solver. The discretized species, turbulence and aero¬ 
sol moment equations result in a block tridiagonal coeffi¬ 
cient matrix which is solved using a block tridiagonal 
solver. All four sets of difference equations are solved in 
delta form to reduce the number of inversions and thus 
the computational effort required to converge non-linear 
iterations. In delta form, the system of discretized non¬ 
linear equations 

[A]{ x y +l = (by ( 35 ) 

can be written as 

[A]{Axy +l =( r y (36) 

where [A] is a square coefficient matrix, {x}" +1 is the solution 
vector, {Ax}' 2+1 is the change in solution vector in the delta- 
form, n is the iteration level and ( r) n = (b) n — [A]{x} n is the 
residual vector evaluated at the previous iteration. In this 
scheme, the residual goes to zero for the converged solution 
vector {x} for all coefficient matrices [A]. However, since 
both linear and nonlinear terms in the system of governing 
equations are treated implicitly, the elements of the coeffi¬ 
cient matrix [A] depend on the solution, time-step, and grid 
related metrics terms, and thus the matrix needs to be 
inverted occasionally throughout the solution procedure to 
maintain effective convergence. Between matrix inversions, 
a computationally inexpensive back-substitution procedure 
is used using the previous LU decomposition for each system 
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Fig. 3. Cross-plane (/', k) equation locations and groupings for 2-D internal and external and quasi 1-D flows. C = continuity, £ = ^-momentum, r\ = rj- 
momentum, £ = ^-momentum, be = boundary condition. Subscripts denote the unknown associated with each equation in coefficient matrix. Fluid state 
and energy equations as well as all moment and species equations are located at the grid points. 


of equations. Whenever the character of the solution 
changes during the marching procedure (due to a change 
in boundary conditions for instance), the matrix must be 
inverted once again. 

4. Code validation: comparison to limiting case solutions 

A number of numerical simulations were carried out to 
validate the aerosol dynamics mechanisms in the model. 
These are chosen based on available limiting cases found 
in the literature for individual particle dynamics phenome¬ 
non important in strongly coupled gas particle flows. These 
include droplet coagulation, droplet condensation and 
evaporation, and combined nucleation, coagulation, and 
condensation. Finally, a validation case for a fully coupled 
case is presented. Validation of the bulk flow, turbulence, 
and particle transport equations have been carried out in 
our previous work [41,42] and are not revisited here. Rele¬ 
vant properties used in the following calculations are listed 
in Table 2. Grid convergence is here defined as a less than 
1% change in the outflow variable values with a 50% 
increase in the number of grid points. 

4.1. Droplet coagulation in uniform flow 

Comparisons with [55] and [56] for aerosol distribution 
properties for particles undergoing coagulation at finite 
times and as time approaches infinity were performed by 
running pseudo 1-D cases. Finite time conditions were 
achieved by setting a slip boundary condition at the wall, 
in a uniform flow and then setting U= 0.025 m/s for vari¬ 
ous mean particle radii and initial standard deviations of 


the particle size distribution. For uniform flow for 1 m, a 
residence time of 40 s is obtained. Initial particle size distri¬ 
butions of (Jqo = 1 and 2 were chosen with an inlet particle 
number concentration of 1 x 10 15 #/m 3 . For all prescribed 
conditions, cr g approaches the theoretical value of 1.355 for 
free-molecular limit and 1.32 for continuum limit as time 
approached infinity. Cases were run on a number of grids 
and grid convergence was reached with 101 axial grid 
points. Figs. 4 and 5 show comparisons to the results for 
finite time cases with the theoretical predictions. Agree¬ 
ment with the current model is excellent. Results show that 
larger standard deviations cause more rapid particle coag¬ 
ulation and thus more rapid increases in mean particle 
diameter. 

4.2. Water condensation and evaporations in uniform flow 
in the low concentration limit 

Fig. 6 shows comparisons of the current method with 
the theoretical predictions of [47] for particle condensation 
and evaporation in uniform flow. Flow conditions are as in 
the previous coagulation cases. Initial particle concentra¬ 
tions are set to 1 x 10 5 #/m 3 with initial standard devia¬ 
tions of the particle size distribution equal to unity. For 
condensation, the initial particle size is 10 pm with a super¬ 
saturation of water vapor equal to 2.023 and for evapora¬ 
tion, the initial particle size is 20 pm with a supersaturation 
of water vapor equal to 0 (dry air). Cases were run on a 
number of grids and grid convergence was reached with 
61 axial grid points. Agreement between calculated and 
theoretical particle sizes is good, though the model slightly 
underpredicts condensation as compared to the theory. 
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Table 2 


Properties of gases and liquids used in calculations 


Value 

h 2 o 

Air 

CsOH 

Units 

MW 


28.964 

149.91 

g/mol 

7 

1.40 

1.40 

1.40 

ND 

Pr 

1.00 

0.71 

0.71 

ND 

Pp 

1000 

N.A. 

3675 

kg/m 3 

E = A + BT + CT 2 + DT 3 

A = 0.05748778 

NA 

A = 0.15286 



B = 3.02908538 x 10~ 4 


B = —1.0660 x 10 -4 



C= -1.0642578 x 10" 6 


C= 2.190 x 10“ 8 



D = 7.0700729 x 1(T 10 


D = 0.00 


C p = A + B/T 1 + C/T+ DT+ 

A = 19.4174489 

A = 25.84546229 

A = 50.92954195 

J/mol/K 

ET 2 + FT 3 + G/7 4 + H/T 5 

B = —228110.480 

B = 47860.'76772 

B = —223544.1081 



C= 3110.1366 

C=0.00 

C=0.00 



£> = 0.021859236 

D = 0.008900406 

D = 0.004977335 



E = —1.522801 x 1(T 06 

E= -2.34795 x 10“ 6 

E=- 5.82599 xl0“ 7 



F= -1.73909 xlO" 09 

F= 2.1164 x 10 -10 

F= 0.00 



G = 4.9167 x 1(T 13 

G = 0.00 

G = 0.00 



H= -3.89863 x 10" 17 

H = 0.00 

H = 0.00 


Ps 

610.483 ( 27 y' 2 ) 51409 x Exp(24.974 T ~ 22i - 2 ) 

NA 

101,325 x Exp(2.3025809 + 4.620)) 

Pa 

(7 

2.641 

3.711 

5.06 

m- 10 

Zk 

809.1 

78.6 

1129 

K 



Location (m) 

Fig. 4. Comparison of particle evolution in uniform flow due to coagulation in the free-molecular regime with the model of [55]. 


4.3. CsOH nucleation, condensation and coagulation in the 
low concentration limit 

To validate the particle nucleation and condensation 
terms and gas phase species coupling, comparisons are 
made with the sectional model of [23] in the one-way cou¬ 
pling limit, i.e. the particles are assumed not to affect the 


gas properties, though the condensing species was con¬ 
served. Thus the fluid field was solved and the gas phase 
species and aerosols where solved using these frozen veloc¬ 
ity and temperature fields. In this case, a mixture of CsOH 
vapor and air is introduced into a one-dimensional temper¬ 
ature field. Calculations where performed using the pre¬ 
scribed temperature field. All other parameters were 
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Location (m) 

Fig. 5. Comparison of particle evolution in uniform flow due to coagulation in the free-molecular regime with the model of [56]. 



Location (m) 

Fig. 6. Comparison of condensation and evaporation rates with the model of [47]. 


allowed to vary. Inlet conditions are T = 520 °C (793 K) 
ramping to 480 °C (753 K) at x=.lm, U = 1 m/s and 
P= 1 atm. Inlet CsOH concentration is 1 g/N m 3 . This 
results in a CsOH mass fraction of 8.4 x 10 -4 . The correc¬ 
tions of [34,17] where used in the calculation of the conden¬ 
sation and nucleation rates along with the properties in 
Table 3. The resulting number concentration of particles 
is low enough that coagulation is unimportant and so 
nucleation and condensation dynamics can be isolated. 


Unlike the previous cases, due to the very low flow rate 
and high nucleation rate, gradients in particle properties 
are extremely sharp and the aerosol equations were found 
to require finer grids due to the high nucleation rates and 
resulting high gradients in the aerosol moments and a large 
number of grid points were needed to achieve grid conver¬ 
gence. This case was computed using a number of equally 
spaced grids ranging from 201 axial points. Grid conver¬ 
gence for the particle dynamic equations was reached at 
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Table 3 


Comparison of the results for CsOH particle evolution in a temperature 
gradient to the model of [23] 



Current lognormal 
model (3 eqns) 

Discrete sectional 
[23] (300 eqns) 

dp g (pm) 

0.147 

0.147 

M 0 (#/N m 3 ) 

1.2 x 10 14 

1.2 x 10 14 


1.11 

1.22 

^csoh (Peak) 

3.25 

3.26 

% Condensed CsOH 

71.9 

71.8 


3001 points. Resulting fluid and particle properties are 
shown in Fig. 7 as a function of axial location. The satura¬ 
tion ratio reaches a peak value of approximately 3.26, at 
which point large scale nucleation occurs and CsOH is rap¬ 
idly converted from gas to condensed phase. The outflow 
value of the gas phase mass fraction was 2.36 x 10 _4 kg/ 
kg resulting in 71.9% conversion for CsOH vapor to the 
condensed phase. Comparisons with the sectional model 
of [23] using 300 particle size bins for the same conditions 
are shown in Table 3. The lognormal model and the sec¬ 
tional model were found to agree quite well as a signifi¬ 
cantly reduced computational cost, though the sectional 
method predicts a wider particle size distribution. 

4.4. Condensation shock formation in a two-dimensional 
wet steam nozzle 

Validation of coupled flow and aerosol behavior 
requires simultaneous comparison of gas phase and aerosol 
phase parameters. Experiments that report both aerosol 


and flow characteristics concurrently are scarce in the liter¬ 
ature. Comparisons are made with the experimental results 
of [57] for condensing water aerosols in a de Laval nozzle 
using wet steam. In that work, pressure and d 32 (Sauder 
mean) aerosol size were simultaneously measured along 
the centerline of the nozzle. Water vapor enters the nozzle 
subsonically and at a pressure below the saturation pres¬ 
sure. As the flow accelerates, pressure and temperature 
decrease, creating supersaturation conditions in a region 
downstream of the nozzle throat. Particles first form by 
homogeneous nucleation and then grow by condensation 
and coagulation. As mass is transferred from the gas phase 
to the liquid phase, heat is released forcing a “condensation 
shock” in the diverging section. The strength and location 
of the shock is directly related to the rates of mass and heat 
transfer due to nucleation and condensation, as well as 
indirectly due to changes in the aerosol properties due to 
the mass loading of droplets in the flow. These influences 
then couple back to the bulk flow which in turn alter the 
local thermodynamic conditions, the growth of the bound¬ 
ary layer, the gas expansion rate and the particle nucleation 
and growth rates. 

The aerosol governing equations are solved in one- 
dimension and two-dimensions (inviscid and viscous/tur- 
bulent) under the conditions described in [57]. 1-D 
simulations were performed to assure £-grid independence. 
The axial grid was then used for 2-D calculations and 
r]-g rid independence studies were run. Grid converged solu¬ 
tions were achieved on a 175 x 45 grid using symmetry con¬ 
ditions at the centerline. For these cases, the total pressure 
at the inlet is 25 Torr. The inlet Mach number corresponds 



Fig. 7. Fluid, vapor and particle properties for CsOH droplet nucleation in a temperature gradient. 
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Fig. 8. Comparisons between the present 2-D viscous/turbulent, 2-D in viscid, and 1-D results with the experimental results of [57]. 2-D viscous/turbulent 
(P/Pt =_i_; dg = -A--), 2-D inviscid (P/Pt = _2. ; dg = -B-), 1-D (P/Pt = .3.; dg = -C-). A' Nuc = 2.5 x 10 3 , K Con = 0.1, K z = 1.25. 


Variation in Fluid and Aerosol Properties along Centerline 
Nozzle C of Moore et 5 (1973) 



Fig. 9. Variation in non-dimensional fluid (Mach number ( M ), total 
enthalpy ( H t ), temperature (F f ) and pressure (P)) and aerosol properties 
(transformed aerosol number (Wo) and volume (W\) concentrations) 
along the center line for [57] nozzle C (viscous/turbulent). 


to choked flow at the throat based on isentropic conditions. 
This is valid since significant aerosol formation does not 
occur under the prescribed conditions until after the throat 
where Mach number is greater than unity and no signifi¬ 
cant upstream influence exists. The inlet temperature corre¬ 
sponds to 21 K below the temperature required for 
saturated dry steam for these conditions (P* f = 359 K). 
The Reynolds number, based on half channel height and 
inlet conditions, is 2.06 x 10 5 . The resulting flow is turbu¬ 
lent. The condensation and nucleation rates and surface 
tension are adjusted based on the viscous/turbulent solu¬ 
tion (Af Nuc = 2.5 x 10 3 , K Co n = 0.1, and K x = 1.25). 

As shown in Fig. 8, the results of the present model 
agree with data quite well for both gas phase and particle 
phase quantities. The condensation shock due to aerosol 
formation is captured in both strength and location, 
however, [57] only reported the d 32 aerosol size at one 
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Fig. 10. Comparison of mixture Mach number (Ma) for viscous/turbulent and inviscid cases for [57] nozzle C. 
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Fig. 11. Comparison of non-dimensional fluid temperature (7» for viscous/turbulent and inviscid case for Moore et al. (3) nozzle C. 
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Fig. 13. Comparison of aerosol volume concentration (Ml) for viscous/turbulent and inviscid cases for Moore et al. (3) nozzle C. 


position. The case is also modeled using an inviscid wall tities) and in quasi-l-D. It can be seen that the effects of 
(with slip boundary conditions on flow and aerosol quan- viscosity are very strong. Fig. 9 shows the centerline flow 
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and aerosol properties for the turbulent calculation. As 
particles form in the nozzle, energy is released into the 
gas phase causing enthalpy, temperature and pressure to 
rise. The aerosol formation process is eventually halted 
as the dropping pressure (due to expansion and scaveng¬ 
ing of the gas phase) overcomes the effects of temperature, 
thus decreasing the saturation ratio below the critical 
level. However, with additional surface area available in 
the particle phase, a significant amount of mass and 
energy transfer continues to take place as particles grow 
by condensation. The viscous boundary layer (Fig. 10) 
slows the flow, and causes the temperature to rise to stag¬ 
nation conditions at the wall (Fig. 11). The elevated tem¬ 
perature in the vicinity of the surface locally reduces the 
saturation ratio and suppresses aerosol formation 
(Fig. 12). This reduces the rate of heat and mass transfer 
between phases and changes the location of the condensa¬ 
tion shock. The resulting aerosol volume concentration is 
shown in Fig. 13. Note the resulting low aerosol concen¬ 
tration near the wall for the viscous case as opposed to 
the high concentration in the inviscid case. 

5. Conclusions 


flow properties can be strong and multi-dimensional effects 
are found to be particularly important. Computational 
experiments reveal that the influence of viscous boundary 
layers are vital when predicting heat and mass transfer 
effects even in relatively high Reynolds number turbulent 
nozzle applications. This is due to the exponential depen¬ 
dency of the nucleation rate on the saturation ratio. This 
influence is found to be prominent even beyond the near 
wall viscous boundary layer region. The damping of nucle¬ 
ation significantly reduces the availability of surface area 
for heat and mass transfer by condensation and changes 
the overall expansion, nucleation and coagulation rates in 
the core of the flow as well as the strength and location 
of the condensation shock. Overall, these results indicate 
that an efficient moment based Eulerian approach com¬ 
bined with an economical RNS solution procedure allows 
the simultaneous resolution of important details of the 
bulk flow, the polydisperse particle properties and the 
interphase coupling effects in strongly coupled particle 
laden flows. 

Appendix A. Non-Dimensionalization 


A moment model for polydisperse aerosol behavior 
incorporating the effects of nucleation, coagulation, evapo¬ 
ration, condensation, convection, diffusion, inertia and 
thermophoresis has been developed in general non-orthog- 
onal coordinates and with the ability to capture strong 
mass and energy coupling between the gas and condensed 
phases in multi-dimensional flows. Limiting cases for parti¬ 
cle evaporation and condensation and coagulation in the 
free molecular and continuum regimes have been modeled 
and validated against solutions in the literature. Combined 
nucleation/condensation dynamics have been modeled and 
validated against a sectional aerosol model for the dynam¬ 
ics in the low concentration limit where fluid particle cou¬ 
pling is one-way. The validated model has been successfully 
applied to water droplet formation in converging diverging 
nozzles with strong heat and mass transfer between the gas 
and particle phases and has been shown to be able to sim¬ 
ulate both the polydispersity of the particle size distribu¬ 
tion and mass and energy coupling effects between 
phases. Comparisons with experimental results for conden¬ 
sation shock formation in a two-dimensional wet steam 
nozzle show good agreement both in terms of the size of 
the droplets and the strength and location of the condensa¬ 
tion shock. It has been demonstrated that the present 
moment model is useful in solving for the evolution of 
polydisperse aerosols in general flow applications. Aerosol 
properties such as the geometric mean particle size, mass 
and number concentrations, and the standard deviation 
of the particle size distribution can be efficiently calculated 
directly without having to sum up results for monodisperse 
aerosols over the entire particle spectrum and with a min¬ 
imum computational effort. As these results indicate, the 
two-way coupling effects between aerosol behavior and 
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Appendix B. Properties of the lognormal distribution 
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